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Abstract. - The nonlinear dynamics of an elastic filament that is forced to rotate at its base is 
studied by hydrodynamic simulation techniques; coupling between stretch, bend, twist elasticity 
and thermal fluctuations is included. The twirling-overwhirling transition is located and found 
to be strongly discontinuous. For finite bend and twist persistence length, thermal fluctuations 
lower the threshold rotational frequency, for infinite persistence length the threshold agrees 
with previous analytical predictions. 



The dynamics and morphology of elastic filaments are of interest in various fields en- 
compassing systems on many different length scales [1]. Relevant examples are provided by 
biopolymes and bio-assemblies, including DNA, actin filaments, microtubules, and multicel- 
lular organisms like Bacillus subtilis [2]. Modern micro-manipulation techniques allow to 
observe and analyze filament dynamics on the single- molecule level [3], prompting a detailed 
understanding of the combined effects of fluctuations, hydrodynamics and elasticity. Static 
properties of elastic rods under external force loads have a long history of study, dating back to 
Euler [4] . Of current interest is the time-dependent behavior of such driven filaments and the 
possibility of shape bifurcations, especially on the biologically relevant nano-to-micro length 
scales, where viscous hydrodynamic dissipation dominates over inertia. 

In this Letter, we describe a hydrodynamic simulation technique for elastic filaments with 
arbitrary shape and rigidity and subject to external forces or boundary conditions, includ- 
ing full coupling between thermal, elastic and hydrodynamic forces in low- Reynolds- number 
flow. We apply our method to the whirling dynamics of a slender filament that is axially 
rotated at one end at frequency uo (while the other end is free). This model system has 
first been studied analytically by Wolgemuth et al. [5], who showed that a critical frequency 
uj c separates whirling (steady-state crankshafting motion with axial spinning) from twirling 
(diffusion-dominated simple axial rotation) by a supercritical Hopf bifurcation (i.e., a con- 
tinuous shape transition). Subsequently, a zero-temperature simulation has been performed 
using the immersed boundary method [6], where the microscopic structure of the filament 
was modelled by interconnected springs. In contrast to analytic predictions, the filament 
was shown to undergo a subcritical (i.e. discontinuous) shape transition from twirling to a 
strongly bent state where the filament almost folds back on itself (termed ov er whirling) . The 
transition frequency uo c was found to be smaller by a factor of 3.5 compared to the analytic es- 
timate. Using our simulation technique, we first map out the stability diagram in the absence 
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of thermal fluctuations as a function of the driving rotational frequency; u c is identified as the 
instability point of twirling and agrees quantitatively with analytical prediction [5]; however, 
a stable small-amplitude whirling state is absent and instead a strongly bent overwhirling 
state is found for all frequencies uo larger than oj c . In the presence of thermal fluctuations the 
transition frequency cu* decreases as the filament becomes more flexible and the discontinuous 
nature of the transitions is manifested by pronounced hysteresis. 

To proceed, consider a filament with circular cross section and contour length L, parame- 
terized by arclength s. A generalized Frenet orthonormal basis {£1,62,63} is defined at each 
point of the filament centerline r(s), where 63 points along the tangent and 61, 62 lie in the 
plane normal to 63 such that the basis forms a right-handed triad, i.e. 62 = 63 x 61. The 
strain vector ft(s) = (^1,^2^3) characterizes the shape of the filament through the kine- 
matic relation d s ej = fi x 6j, where n — (ftf + ft^) 1 ^ 2 * s the curvature and ^3 the twist 
density. The elastic energy of an elastic inextensible rod reads [4] 

E =j£ dsKAQ,) 2 + (Afi 2 ) 2 ] + ds(AV 3 ) 2 , (1) 

where A and C are bend and twist rigidity constants, and Aft = fl — ft . The ground state 
shape of the filament is thus controlled by the intrinsic curvatures f2j, ^2 an d the intrinsic 
twist ^3. According to the principle of virtual work, the force per length acting on the 
filament is obtained via f = —SE/Sr keeping the rotation angle about the tangent zero, i.e., 
Scj) = 62 • <56i =0 [7]. Likewise the torque about the tangent is obtained via m = —SE/Scf) 
without moving the centerline, i.e., 5r = 0. In the low- Reynolds- number limit, the local 
force and torque balance with the viscous drag (long-range hydrodynamic interaction will be 
considered later), leading to the equations r = /if and <j) = /i r m, where r = dv/dt and /i, \i r 
are the translational and rotational mobility constants. 

For an efficient simulation code, it is crucial to introduce a laboratory- fixed reference frame 
{x, y,z} and to parameterize ft (thereby the energy E) in terms of r [8]. The unit vector 
n = txz/|txz|isby construction normal to the tangent t = 63. Taking the binormal b = 
t x n, we find a local right-handed triad {n, b, t} at the filament centerline position r(s). Since 
(61, 62) and (n, b) are coplanar by definition, they are related as 61+262 = exp[— z#(s)](n+ib). 
The rotation angle about the tangent thus satisfies 5<j) = 56 + b • <$n, indicating that r and 
are not independent in computing the translational force f . Using the kinematic relation 
Qi = — 62 • £^63, ^2 = 61 • £^63, and Q3 = 62 • <9 s 6i, it is straightforward to express ([T|) in 
terms of r and 0. 

In our dynamic simulation, the filament is modelled as a chain of N + 1 connected spheres 
of diameter a. Each bead is specified by its position rj and the twist angle with respect 
to the laboratory frame, 0j, from which the strain vector flj is calculated at each sphere 
point. The unit tangent is thus given by tj = Uj/|uj|, where Uj = r J+ i — rj is the bond 
vector. The simplest symmetric discretization of ^3 is by = [^( e 2 * d s ei — e± • d s e2)]j = 
fj sin(#j — Oj-i) + gj cos(0j — #j-i), where / and g are functions of {r^} only. The total elastic 
energy E involves the stretching contribution: E tot = E[r, 6] + XljLV ^/(^ a )(\ u j \ ~ °) 2 - For 
an isotropic rod, the stretching modulus K is given by K = 16 A/ a 2 . Writing the angular 
velocity as dt<fi = d t + b • d t n and using V Yj = — b • V rj n, we arrive at the coupled Langevin 
equations 

AT+l 



iV+l 



^v 3 E to t\ Q + ^o k E t ot\ Y ^k • V rj n fc 



k=l 



(2) 
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0i(t) = -hi ■ hi - HrVe 3 E to t + Si(t). (3) 

Hydro dynamic interactions between two spheres i and j are included via the Rotne-Prager 
mobility tensor /Xij(r^) = 1/(8^777^ )[1 + + a 2 /(2r^)(l/3 - r^r^)], where = r^/r^ 
and 77 the solvent viscosity [9]. For the translational and rotational self- mobilities of the 
spherical monomers with diameter a we use fin = l/(37rr/a) = fiol and fi r « l/irrja 3 [10]. 
The hydrodynamic coupling of rotation and translation between two spheres, the so-called 
ro^/et effect [10], decays fast in space (~ 1/r 2 ) and is neglected, similar to previous studies [5]. 
Rotlet corrections are examined separately [11]. The vectorial random displacements £(t) 
and S(t) model the coupling to a heat bath and obey the fluctuation-dissipation relations 
= 2k B Tn ij 5(t-t), {EiifiZtf)) = 2k B T^ r 5 ij 5{t-t') and (~^)^')> = 0, 
which are numerically implemented by a Cholesky factorization [9]. 

For the numerical integrations we discretize Eqs. ([2]) and (j3]) with a time step A and rescale 
all lengths, times and energies, leading to the dimensionless parameters A = AksTfio/a 2 , 
/ /Jo — 3, L p /L = A/ \aNksT) and u = coa /(/iofceT). We set the twist-bend 
rigidity ratio to C/A = 2. This choice, which corresponds to a negative Poisson ratio a = 
A/C — 1 [4], is relevant for typical biopolymers [12, 13]. We also studied the C/A = 1 case, 
with no qualitative difference [11]. For sufficient numerical accuracy we choose time steps in 
the range A = 10 -4 -10 -7 . Output values are calculated every 10 3 -10 4 steps, total simulation 
times are 10 6 -10 8 steps. Clamped boundary conditions at the forced end, i.e. <9 s r(0) = 0, are 
realized by fixing the first two monomers in space by applying virtual forces, which also act 
(via the mobility tensor) on the rest of the filament. A mobile filament is obtained by switching 
off the virtual forces along the axis (chosen as the x direction). The rotational driving at the 
base imposes 0\(t) = uj. Force- and torque- free boundary conditions are adopted for the other 
end. As initial shape of the filament we take circular sections specified by the angle a; see the 
inset of fig. [T] (g). The straight rod configuration thus corresponds to a = 0. The number of 
beads is in the range L/a = N = 20 — 40. 

At low rotational frequency, uo < uo Cl the rod is twisted but remains straight and the torque 
at the base balances the total rotational drag, \i~ x ujL ~ CQs(0). On the scaling level, the 
rod buckles when the twisting torque C£^3(0) becomes comparable to the bending torque, 
A/L, giving a critical frequency uo c ~ fi r A/L 2 independent of the twist rigidity C. The 
asymptotically exact linear analysis based on slender-body hydrodynamics predicts a critical 
frequency [5] 

co c ^ 8.9/i r A/L 2 = S.9fi r k B T(L p /L) 2 /L p . (4) 

The time evolution of the rod shape for TV = 30, persistence length L p /L = 10 3 for uo/uo c = 1.20 
(starting with a straight shape a = 0) is shown in fig. [T^i. Twirling is unstable against thermal 
disturbance and the filament buckles to relieve twist, leading to crankshafting motion. As will 
be discussed in more detail below, the radial distance of the free end from the rotational axis, 
called i?, initially increases exponentially in time until the filament bends over (see fig. [TJd) 
and the steady overwhirling shape is reached. Overwhirling is a combination of rigid-body- 
like rotation with frequency x an d axial spinning. The crankshafting frequency x increases 
drastically across the kinetic transformation state, while the torque applied at the base, M, 
shows a small but steep decrease (see fig. Cfc-d); in other words, the transition to the whirling 
state can be viewed as a way to reduce the dissipated power. In figs. [1^-1 and e-2 the steady 
state values of x an d M are plotted against uo/uo c across the transition for a rather stiff rod 
with L p /L = 10 3 . A pronounced hysteresis is revealed, suggesting that the observed transition 
is strongly discontinuous and the final state depends not only on the driving frequency uo but 
also on the initial shape of the filament (described by a). This stands in vivid contrast to the 
analytics predicting a continuous transformation [5]. The bifurcation frequencies depend on 



4 



EUROPHYSICS LETTERS 




Fig. 1 - (a) Sequence of snapshots of a rotating filament of length N — 30 and stiffness L p /L = 10 3 
for uo/uoc = 1.20. Corresponding time evolutions of (b) the end-point radius R, (c) the crankshafting 
frequency % and (d) the torque applied at the driving end M. (e-1) Hysteresis loop of x f° r a iV = 30 
and Lp/L = 10 3 rod. (e-2) same for M. (f-1) Observed transition frequency to* for varying stiffness 
Lp/L for N = 30 rod. (f-2 and f-3) Corresponding time series of x an d w close to the transition 
for Lp/L = 10. (g) Zero-temperature stability diagram as a function of initial bending angle a and 
driving frequency uj for a N — 30 rod. All data are for an immobile filament. 



the rate with which the frequency is changed, but only for very flexible rods does the hysteresis 
disappear. We construct the stability diagram shown in fig. [T] (g) in the (lj, a) plane at zero 
temperature (i.e. the deterministic noise- less limit where L p /L = oo) by following the time 
evolution of a rod with an initial circular shape with bending angle a and the stationary 
twist profile. Twirling is the final stationary state for a filament with an initial shape bending 
angle below the separatrics, while above the separatrics the overwhirling state is obtained. 
The separatrics terminates for a = at uo = uo c as given by eq. (|4j) [5], indicating that uo c 
is an instability point and actually the correct critical frequency at zero-temperature above 



H. WADA and R. R. NETZ: NON-EQUILIBRIUM HYDRODYNAMICS OF A ROTATING FILAMENT 5 

which the filament is unstable against infinitesimal disturbance. By construction, the stability 
diagram does not depend on the stiffness of the filament. The cause of the discrepancy to the 
previous simulation [6] at zero temperatures which yielded a much lower transition frequency 
is unclear. 

At finite temperature (i.e. finite L p /L), the overwhirling transition is observed below uo c : 
In fig. [U (f-1) the transition frequency uo* is plotted as a function of L p /L, obtained from 
simulations where the driving frequency uo is increased in steps of size Auo = uo c / 100 every 10 7 
simulation time steps. The time-dependent behavior of x and uo across the transition is plotted 
in fig. [1] (f-2) and (f-3), respectively, for L p /L = 10, from which uo* is determined. An even 
slower frequency increase rate gives a slightly lower u;*, so the results in fig. [1] (f-1) correspond 
to an upper limit of the transition in the stationary limit. A simple scaling arguments explains 
the observed trend: Balancing the bending energy E ~ Aa 2 / (2L) with /c#T/2, one obtains a 
spontaneous mean squared bending angle (a 2 ) ~ L/L p , which is larger for smaller L p and thus 
leads to reduced stability (see the stability diagram). Below L p /L ~ 5, the filament switches 
frequently between the two states, and the transition is washed out. The dependence of the 
twirling-overwhirling transition on the persistence length could be experimentally tested with 
actin or tubulin filaments of varying length. 

For a stiff rod and short times, we expect universal shape dynamics for frequencies slightly 
above uo c . Consider a small amplitude whirling motion (R/L <C 1) at crankshafting frequency 
X for supercritical frequency e = (uo — uj c )/uj c > 0. As we will show, the small amplitude 
whirling is unstable for any e, but lasts for a long time for small enough e due to critical 
slowing down. The rod shape is approximated as r(s) ~ sx + r^(s), and r(s) = x x x r ±( 5 ) f° r 
crankshafting. Since x grows much more slowly in time than the end-point radius i?, it can be 
estimated by assuming that the transverse drag force per length, xkj-l/^o? is roughly equal to 
the bending force per length, A\r±\/L 4 , giving x ~ Afio/L 4 ~ (a/2L) 2 uo c . After a short time 
transient, the twist density ^3 has built up the stationary linear profile, ^3 = uo/fi r C(s — L), 
and thereby reached ^3 = 0. The difference in rotational velocities about the local tangents 
at s = and s = L, Aa;, satisfies —Auo = x [1 — x • t(L)] (for a derivation see ref. [5]). 
Considering the local balances of viscous and elastic twisting torques \i~ x uoL ~ CQ3 at s = 
and s — L, we see that the net torque AM ~ AuoL must contribute to the whirling 
motion of the filament. On the other hand, the elastic energy E of the rod with the amplitude 
R is estimated to be E ~ 2AR 2 /L 3 , while the power dissipation Pd by the rigid-body like 
rotation is Pd ~ v 2 L/Z\±§ = x 2 ^ 2 £/3/io, where v = xR 18 the rotational velocity. The 
power balance condition, according to which the sum of the change of the elastic energy 
per unit time and the hydrodynamic rigid-body dissipation balances the net torque injected 
into the rod per unit time, implies E + Pd ~ uoAM. Note that the remaining power input 
uj(M — AM) is independently balanced with the dissipation by axial rotation. Approximating 
1 — x • t(L) w 2R 2 /L 2 we finally arrive at R(t) ~ exR(t). The end point radius R(t) thus 
evolves exponentially in time for e > 0, and the growth rate T(uj) = d(\ogR)/dt satisfies the 
linear relation T/x ~ (uo — uo c ) /uo c . In fig. [2] (a), the numerically obtained T/x, which is almost 
constant while R/L <C 1, is plotted as a function of the super cr it icality e = (uo — uo c )/uo c for 
filaments of various length and stiffness. The thermal noise is switched off for improved data 
accuracy. The data scale linearly for small e in agreement with the prediction (with a numerical 
prefactor of the order of unity), verifying the validity of the arguments and approximations 
presented above. For larger £, nonlinearity comes in, leading to deviations from the simple 
linear law. For the numerical calculations the presence of this critical slowing down can be a 
problem because stationary states are only obtained after a long waiting time. 

In the case when the filament is allowed to move along the axis of rotation, a finite average 
velocity is observed in the overwhirling regime: the filament is propelled due to the formation 
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Fig. 2 - (a) Growth rate of end radius F divided by the crankshafting frequency % as a function of 
the supercriticality parameter e = (00 — oo c )/oo c for rods with various length and stiffness obtained in 
the noise- less limit, (b) Rescaled propulsion velocity along the rotation axis V\\ = V\\a/ fioksT in the 
overwhirling state, plotted against co/oo c (Inset: actual time evolution of the base position x\ across 
the shape transition for e = 0.1). (c) Vy versus the rescaled external force F ext = aF ext /kBT for 
e — 0.1. (d) Efficiency 77 as a function of F ex t with the parabolic fit (dashed line). All data in (b)-(d) 
are for N = 30 and L p /L = 10 3 (i.e. at finite temperature). 



of a helical shape which breaks the time-reversal symmetry of the problem. Accordingly, a 
jump in the propulsion velocity along the rotation axis V\\ is observed across the overwhirling 
transition, resulting in a drastically amplified forward thrust of the filament (see the inset 
of fig. [2] (b)). The propulsion velocity V\\ in the overwhirling regime is plotted in fig. [2] 
(b) as a function of uj/uj c . Note that the overwhirling state is metastable for a wide range of 
frequencies below cj* and thus the velocity is non-zero. The non-monotonic increase of V\\ with 
ou implies a complicated shape change of the rod as a function of a;. To study the swimming 
efficiency of this self-propelling object, we apply an external force F ext at the filament base. 
The power converter efficiency is the ratio of the propulsive power output and the rotary 
power input, 77 = F ex tV\\/{Muj), where a positive F ext is defined such as to work against 
the natural propulsion direction of the filament. The rotational and translational mobilities 
of the whole filament are defined by uj = fi rr M + {i r tF ex t and V\\ = Li tr M + UttFext [14]. 
Although for perfectly stiff propellers the mobility matrix is constant and symmetrical (i.e., 
in rt = /i tr ), due to the flexibility of the filament, the symmetry here is broken and the mobilities 
depend on the external torque M and force F ext in a complicated non-linear manner. In the 
simulation, a substantial forward thrust is observed independent of the rotation direction, 
which suggests that \i tr changes its sign when the torque N is reversed, similar to results 
obtained for a rod with vanishing twist rigidity that rotates on the surface of a cone [15]. V\\ 
varies linearly with F ext as shown in fig. [2] (c), which means that ji u is independent of the 
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external force F ext . Neglecting /i r t, which is actually vanishingly small, while keeping /i tr , we 
obtain r] = (n rr n tt Fext + Mtr^^ct)/^ 2 , indicating that the efficiency becomes parabolic as a 
function of the external force. Figure [2] (d) shows the numerical data of r] with the parabolic 
fit for e = 0.1, which describes the date quite nicely. The highest efficiency is only of the order 
of 0.01%, which means that this self-propelling filament is quite inefficient compared to other 
known examples involving helices and other chiral objects [15], since most of the input power 
is dissipated by the axial spinning. 

In summary, the elastohydrodynamics of a spinning filament in viscous fluid is studied by 
simulations including full coupling between elastic, thermal and hydrodynamic effects. Quan- 
titative agreement with the analytically predicted critical frequency is obtained [5] in the 
absence of thermal fluctuations, in contrast to numerical simulations of the same problem [6]. 
We give evidence for the discontinuous nature of the twirling-overwhirling transition, in qual- 
itative agreement with previous numerical works at zero temperature [6] but in contradiction 
to the analytical theory which might have to do with the neglect of non-linear effects [5]. 
Thermal fluctuations play a significant role in the transition behavior of a filament and lead 
to a decrease of the critical frequency in a range of filament stiffnesses L p /L that is relevant 
to biopolymers. The parametrization proposed in this paper is entirely based upon a local- 
form description and advantageous for numerical implementations. Our dynamic simulation 
method is easily generalizable to other geometries such as rings or helices. A number of in- 
triguing problems have to do with pulling on flexible nanospring [16] or buckling of twisted 
rods [1,12]. Of particular biological interest is a spinning helical filament. Taking, for exam- 
ple, a diameter value of a « 20 nm and a typical stiffness L p /L « 10 3 of a bacterial flagellum 
at physiological conditions [17], we find uo c ~ 10 3 s _1 in water with the viscosity r] ~ 10 -3 
Pas, which is easily achievable in laboratory experiments. Details on the interplay of thermal 
fluctuations and shape instabilities will be reported elsewhere [11]. 

We thank M. Manghi for valuable discussions and the program for Research Abroad of 
the Japan Society for the Promotion of Science (JSPS) and the German Science Foundation 
(DFG, SPP1164) for financial support. 
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